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Abstract 

We  study  a  diffusive  energy -balance  climate  model, 
governed  by  a  nonlinear  parabolic  partial  differential 
equation.   Three  positive  steady-state  solutions  of  this 
equation  are  found;  they  correspond  to  three  possible 
climates  of  our  planet:   an  interglacial  (nearly  identical 
to  the  present  climate),  a  glacial,  and  a  completely  ice- 
covered  earth.   We  consider  also  models  similar  to  the  main 
one  studied,  and  determine  the  number  of  their  steady  states. 
All  the  models  have  albedo  continuously  varying  with  latitude 
and  temperature,  and  entirely  diffusive  horizontal  heat 
transfer.   The  diffusion  is  taken  to  be  nonlinear  as  well  as 
linear. 

We  investigate  the  stability  under  small  perturbations 
of  the  main  model's  climates.   A  stability  criterion  is 
derived,  and  its  application  shows  that  the  "present  climate" 
and  the  ''deep  freeze"  are  stable,  whereas  the  model's  glacial 
is  unstable.   A  variational  principle  is  introduced  to  confirm 
the  results  of  this  stability  analysis. 

We  examine  the  dependence  of  the  number  of  steady  states 
and  of  their  stability  on  the  average  solar  radiation.  The 
main  result  is  that  for  a  sufficient  decrease  in  solar  radia- 
tion (about  2  percent)  the  glacial  and  interglacial  solutions 
disappear,  leaving  the  ice-covered  earth  as  the  only  possible 
climate . 
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1.   Introduction 

The  concept  of  a  climate  is  one  of  those  abstractions 
which  appears  to  be  self-evident  to  the  layman,  but  is  by 
no  means  well  defined  scientifically.   The  intuitive  idea 
of  a  climate  has  two  aspects: 

(a)  the  most  important  features  of  atmospheric  phenomena; 

(b)  the  average  behavior  of  these  phenomena  over  a  suitably 
long  time  interval  and  over  sufficiently  large  areas. 
The  difficulties  start  when  one  tries  to  give  a  precise 

meaning  to  the  key  words  "most  important",  "suitably  long" 
and  "suitably  large".   We  start  with  "suitably  long" ; 
clearly,  a  year  is  an  absolute  lower  bound  for  a  reasonable 
averaging  time  interval,  since  daily  and  seasonal  variations 
should  be  excluded.   To  decide  over  how  much  longer  than  a 
year  the  averaging  should  be  performed,  one  has  to  look  at 
the  record.   There  are  three  kinds  of  records:   instrumental, 
the  length  of  which  is  of  the  order  of  hundreds  of  years, 
historical,  of  the  order  of  thousands  of  years,  and  geological, 
of  the  order  of  hundreds  of  thousands  of  years.   These 
records  show  that  features  of  the  atmosphere  change  on  all 
the  time  scales  represented  in  them  (e. g. ,  Robinson ,  1971). 
Thus  it  would  appear  at  first  that  it  is  not  possible  to 
distinguish  between  "fast"  variations  in  yearly  averages, 
which  should  be  averaged  out  when  defining  a  climate,  and 
"slow"  variations,  which  should  be  considered  as  "changes 
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of  climate".   Still,  the  geological  record  seems  to  indicate 
that  the  transitions  between  considerably  colder  periods  (ice 
ages  or  glacials)  and  warmer  periods  ("normal"  climates  or 
interglacials)  occurred  over  time  spans  about  ten  times 
shorter  than  the  duration  of  the  relatively  steady  cold  or 
warm  weather  respectively.   This  suggests  what  we  shall  adopt 
here  as  our  operative  definition  of  climate,  viz.,  the  preva- 
lence of  either  warm  weather  (as  we  experience  it  today)  or 
of  cold  weather  (to  mean  a  difference  of  the  order  of  ten 
degrees  centigrade  in  yearly  average  below  the  one  recorded 
in  the  present) . 

We  turn  now  to  the  question  of  which  features  of 
atmospheric  phenomena  are  "most  important" .   Certainly  tempera- 
ture is  one  of  them,  not  only  because  its  changes  left  deep 
traces  in  the  geological  record  (glaciations  in  temperate 
zones,  pluviations  in  the  tropics  —  SMIC,  1971),  but  also 
because  it  affects  all  conditions  of  life  and  because  it  is 
directly  linked  to  the  major  thermodynamic  and  dynamic 
processes  in  the  atmosphere  which  determine  climate  and 
its  changes.   Also,  humidity,  wind  direction   and  intensity, 
cloud  amount,   precipitation,   all   play   a   major 
role  in  determining  what  is  perceived  as  weather  and  hence 
should  be  time-averaged  (and,  possibly,  space-averaged) 
into  climate.   Moreover,  it  is  not  only  the  averages  of 
these  quantities,  but  their  day-to-night  and  season-to-season 
contrast  that  enters  our  intuitive  concept  of  a  climate. 
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Thus,  at  least   their  variance  should  be  included  in  a  more 
complete  mathematical  model  for  climatology. 

In  this  article   we  shall  treat  a  very  simple  model, 
based  on  the  work  of  Sellers  (1969)  and  of  Schneider 
and  Gal-Chen  (1973);  we  hope   that  the  results  will  in 
themselves  be  of  some  significance  for  climate  theory,  as 
well  as  providing  insights  for  devising  and  analyzing  more 
complex  models. 

In  Section  2  the  model  to  be  studied  is  described; 
the  physical  principles  on  which  it  is  based,  as  well  as 
the  empirical  data  it  uses  are  discussed. 

In  Section  3  we  discuss  the  work  of  different  authors 
on  similar  models;  the  similarities  and  differences  between 
their  results  are  pointed  out  and  the  issues  arising  from 
these  results  are  outlined. 

In  Section  4  we  compute  numerically  the  model's 
steady-state  solutions  of  physical  interest,  i.e.,  those 
yielding  positive  absolute  temperatures.   Three  such 
solutions,  corresponding  to  three  distinct  climates  of  our 
planet,  are  obtained:   one  corresponds  to  the  current  climate, 
the  second  to  an  ice  age,  the  third  to  a  completely  ice- 
covered  earth.   In  this  section  vie    also  explore  the  effect 
of  certain  modifications  in  the  model  on  the  number  of 
steady-state  solutions. 

In  Section  5   the  notion  of  stability  for  the  solutions 
obtained  in  Section  4  is  defined  precisely;  it  is  investigated 
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using  a  combination  of  analytical  and  numerical  techniques. 
The  results  are  that  the  present  climate  and  the  ice~covered 
earth  are  stable,  whereas  the  ice  age  of  the  model  is 
unstable . 

In  Section  6  the  effect  of  changes  in  the  solar  radia- 
tion on  the  number  and  stability  of  steady-state  solutions 
is  studied.   The  main  result  is  that  for  a  sufficient 
decrease  in  the  solar  radiation  (about  2%) ,  the  glacial  and 
the  interglacial  solutions  disappear,  leaving  the  ice- 
covered  earth  as  the  only  possible  climate. 
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2 .      The   Model 

The   model    chosen    for   study    is    a    zonally-and-vertically 
averaged  energy -balance    climate   model.      This    means    that 
quantities    in    the   model    are    averaged   over    longitude    and 
height,    leaving   colatitude      $      as    the    only   space   variable. 
The    term  energy   balance   means    that   the   model    is   essentially 
based   on    the   energy   equation    of    fluid   dynamics    and  has    sea- 
level    temperature      u      as    the    only   dependent   variable.       The 
equation   governing   the   model    is 

(la)       C(4')u^   =    R^{(l),u)    -    RQ{(f),u)    +    D(cj),u,u^,u^^)     ; 

C  is  the  heat  capacity  of  the  atmosphere,  land  and  water 
masses;   R.   is  the  heat  absorbed  from  incoming  radiation, 

(lb)  R^   =   Q(<t))  [1  -  a(<}),u)  ]  , 

where      Q   is   high-frequency    solar   radiation    and      a   is    the 
reflectivity    (albedo)    of    the    land   and   sea   surface;      Rq    is 
the   heat   lost   in   outgoing   low-frequency   planetary    re- 
radiation      reaching   outer   space, 

4 
(Ic)  Rq    =    c((t),u)au      ; 

and   D   describes    the    redistribution   of  heat   on    the    surface 
of   the   planet  by   conduction    and   convection, 

(Id)  D   =   ^^  f^[sin    4)    •    k(<^,u)]u^    . 

The  coefficients  and  forcing  terms  in  this  model 
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represent  yearly  averages  of  the  corresponding  quantities  and 
therefore  do  not  depend  explicitly  on  the  time  t.  Averaging 
the  daily  and  seasonal  variations  seems  justified,  since  the 
time  scales  in  which  we  are  interested  in  our  investigation 
are  of  the  order  of  hundreds  and  thousands  of  years.  The  lack 
of  explicit  dependence  on  time  has  the  advantage  that  the 
model  admits,  as  we  shall  see,  steady-state  solutions,  u   =  0, 
which  we  define  to  be  its  climates .   The  purpose  of  this  work 
is  to  study  analytically  and  numerically  the  number  of  these 
climates  and  their  stability  under  perturbations. 

The  first  model  of  this  type,  in  finite-difference  form 
and  without  time  dependence,  was  developed  by  Sellers  (1969). 
The  differential  formulation  is  due  to  Faegre  (1972);  time 
dependence  was  introduced  by  Dwyer  and  Petersen  (19  7  3)  and, 
independently,  by  Schneider  and  Gal-Chen  (19  73) .  Dwyer  and 
Petersen  also  gave  the  outline  of  a  systematic  derivation  of 
(1)  from  the  energy  equation  of  fluid  dynamics,  mentioning  the 
main  assumptions  involved.   An  even  simpler  model  has  been 
proposed  by  Budyko  (1969):  in  it  the  diffusion  term  D  is 
replaced  by  a  nondif ferentiated,  linear  term  in  u,  and  the 
albedo  is  a  simple  step  function  of  u  only;  this  model  was 
also  discussed  very  thoroughly  by  Leith  (1974),  and  by  Held 
and  Suarez  (1974) . 

One  of  the  main  features  of  the  model  (la-d)  is  the  form 
of  the  albedo, 

(2a)        a  =  {b((}))  -  c^  [u^  +  iu-c^z{(^) -uj  _]  }  ^    , 
where  the  meaning  of  the  subscripts  (  ) _  and  {  }   is 
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given    for    a   generic   quantity      h      by 


(2b) 


h      =   min    {h,0} 


and 


(2c) 


0.25,  h    <    0.25, 

h         ,       0.25    <   h    <    0.85, 

^0.85,       0. 85    <    h  , 


The  subscript  c  stands  for  cutoff;  the  cutoff  given  by  (2c) 
eitibodies  the  observed  minimum  and  maximum  values  of  surface 
albedo. 

Snow  and  ice  have  higher  reflectivity  than  bare  ground 
or  water;  since  in  regions  of  lower  yearly  average  tempera- 
ture the  snow  and  ice  cover  persists  for  a  longer  fraction 
of  the  year,  at  lower  temperature  the  yearly  average  albedo 
is  higher;  this  is  expressed  in  the  monotonically  decreasing 
dependence  of  a  on  u.  Further,  the  plausible  assumption  is 
made  that,  above  a  certain  yearly  average  temperature  u  , 
no  snow  or  ice  will  be  present  at  any  time  of  the  year; 

therefore  a  is  independent  of  u  for  u-c^z  >  u   ,  as  seen 

^  2      —     m 

from    (2a)     and    (2b).      The    term      c^z((l))    gives    the   difference 
between    sea-level   temperature    u   and   ground   temperature,    u-C2Z' 

A   serious    drawback   of    the   model    is    that   it   does   not 
include    the   effect   on    the    albedo    of   clouds,    atmospheric 
turbidity,    relative   humidity,    and  vegetation.      The   optical 
properties    of   these    factors    and    their   relationship   to   surface 
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temperature    are    less   well   known    and  cannot  be   easily  para- 
meterized  in   a   model    as    simple    as    the   one   at   hand. 

The  factor   c   in  the  outgoing  infrared  radiation 
term  R-  , 

g 
(2d)  c  =  1  -  m  tanh  (c^u  )  , 

expresses  empirically  the  "greenhouse  effect",  i.e.,  the 

screening  by  the  atmosphere,  in  particular  by  the  clouds 

in  it,  of  infrared  radiation  from  the  earth,  thus  preventing 

part  of  it  from  reaching  outer  space.   Notice  that  c  decreases 

as  u  increases;  this  indicates  that  cloud  formation,  and 

hence  the  opacity  of  the  atmosphere  to  low-frequency  radiation, 

increases  with  increasing  temperature. 

The  function   k(ct),u)  in  (Id)  has  the  form 

c.   -c  /u 
(2e)   k((}),u)  =  k,((j))+k„(<}))g(u)  ,   g(u)  =  -^  e      =  f '  (u)  ; 

-^     "^  u 

k,(d))u^   is  sensible  heat  flux  in  the  atmosphere  and  in  the 
1     (p 

oceans,    whereas      k- (<)>)  g(u)  u  ,       is    latent  heat    flux   in    the 
atmosphere.      Here      k,  ((J)),    kjifti)       are   eddy   dif  fusivities , 
which   parameterize    convective    transports;    true    conduction 
is   known   to  be  negligible    in   the    atmosphere-ocean   system 
on    the   planetary   scale.       In    the   original   Sellers   model, 
k((J),u)      had   the    form 

k((t>,u)    =   kg((J),u)    -    v(u   )  •  (u+f  (u)  )     , 

where  v  is,  for  the  purposes  of  modeling  heat  flux,  mean 
meridional  velocity.   The  theoretical  shortcomings  of  the 
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additional  term  v(u+f(u))   were  pointed  out  by  Robinson 
(1971);  the  practical  difficulties  in  giving  a  good  para- 
meterization are  discussed  by  Sellers  (1973). 

Numerical  studies  of  Schneider  and  Gal-Chen  (1973) 
indicate  that  results  with  the  original  Sellers  model 
(denoted  by  them  as  (S) )   were  very  similar  to  those  with 
the  model  adopted  here  (denoted  in  their  work  as   (SV) ) , 

provided  that  the  numerical  values  of   k  and  k   were 

'^  s 

properly  chosen   (see  the  discussion  on  the  determination 
of  coefficients  further   on).     Furthermore,  the  recent 
work  of  Gal-Chen  and  Schneider  (19  75)  shows  that  there  are 
theoretical  grounds  on  which  the  formulation  (SV)  with  zero 
meridional  velocity  is  to  be  preferred.   These  considerations 
will  be  touched  upon  later,  in  a  different  connection. 

The  constants   c.  ,l£j£5,  u^   ,a,m,as  well  as 
the  empirical  functions  C{<t>),    0(4)),  b(4)),  z{(p),    k^{(p)    and 
k-,((J))  are  made  to  fit  currently  measured  values  of  temperature, 
radiation,  elevation,  albedo  and  heat  flux.   The  functions 
C((t>),  Q((}))  and  z((|))   are  determined  directly  from  measurements. 
The  function  b(({))   and  the  constant   c,  in  (2a)  were 
computed  by  Sellers  (1969)  so  as  to  fit  existing  albedo 
measurements  in  the  northern  and  southern  hemisphere. 
The  constants   m  and  c^  were  also  computed  by  Sellers,  so 
as  to  fit  empirical  data  on  R- ;   a  is  the  Stef an-Boltzmann 
constant.   The  form  of  the  function   g(u)   and  the  constants 
c.,    c_   appearing  in  (2e)  are  derived  from  certain  physical 
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considerations  having  to  do  with  the  thermodynamics  of  wet 

air  and  from  corresponding  empirical  data  (see  Handbook  of 

Meteorology,  1945).   The  functions   k-,((j)),  k^  ((]>)  are  computed 

from  measured  data  on  sensible  and  latent  heat  flux,  k^  {<i))u, 

1     ()) 

and  k„(4))g(u)u,  ,  respectively.  These  computations  are 
based  on  the  measured  temperature  distribution,  denoted 
hereafter  by 

u  ==  u(4))  , 
which  will  be  called  the  data  climate.   Note  that  we  use 
here  the  term  "climate"  only  for  convenience,  instead  of  the 
lengthier  "temperature  distribution";  u(.<^)       is  not  necessarily 
a  steady-state  solution  of  the  model;  we  return  to  this  point 
in  Section  4. 

The  measured  data  are  available  at  intervals  of  10° 
latitude  and  are  given  in  Table  1.   Since  the  previously 
quoted  authors  used  finite-difference  formulations  with  a 
fixed  10°  grid  (except  Faegre  (1972) ,  who  used  a  5°  grid) , 
these  data  were  sufficient  for  their  numerical  work.  In  our 
numerical  work,  however,  variable  grid  size  was  employed, 
and  the  10°  data  were  accordingly  fitted  by  smooth  functions. 
In  fact,  in  order  to  have  an  additional  check  on  the  well 
posedness  of  the  model  (i.e.,  the  continuous  dependence  of 
the  solutions  on  the  data,  commonly  referred  to  in  the 
rreteorological  literature  as  sensitivity)  ,  two  forms  of  curve 
fitting  were  used:   (i)  by  Bernstein  polynomials,  and  (ii)  by 
cubic  splines.   Results  with  the  two  forms  of  curve  fitting 
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were    indeed   very    similar    (see   Table    2)  . 

In    this   work   we    only    consider    syitunetric   solutions    of    (1)  ; 
all   data   are   syitimetrized  with    respect    to   the      equator, 
(|)   =   tt/2.      For   such   data   the    appropriate   boundary   conditions 
are 

(3a)  u^(0)    =    0    ,  (3b)  "^^^Z^)    =    0    ' 

in  the  symmetric  case  these  are  equivalent  to  u  (0)=  u  (11)=  0. 

We  feel  that  the  slight  asymmetry  between  the  northern 
and  southern  hemispheres  could  hardly  have  had  a  major 
influence  on  climatic  change.   Indeed,  the  circulations 
of  the  two  hemispheres  are  practically  separated  from  each 
other  by  the  intertropical  convergence  zone,  which  acts  with 
respect  to  our  model  as  an  insulating  wall.   The  approximation 
involved  in  placing  this  wall  at  the  equator  is  no  worse 
than  other  approximations  in  the  model  (see  also  Held  and 
Suarez,  1974).   A  further  reason  for  symmetrization  will 
become  apparent  in  the  next  section. 
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3.   Previous  Results 

Budyko  (1969)  and  Sellers  (1969)  used  iterative  numeri- 
cal techniques  for  constructing  solutions  of  their  time- 
independent  models.   They  explored  a  range  of  values  of  the 
parameters  appearing  in  the  model,  especially  of  the  solar 
radiation   Q,   and  obtained  one  solution  for  each  set  of 
values .   These  solutions  did  not  depend  smoothly  on  the 
parameter  values;  in  several  instances,  small  changes  in  the 
parameters  led  to  large  changes  in  the  solution.   For 
instance,  an  increase  or  decrease  of  a  few  percent  in  Q 
resulted  in  temperature  changes  leading  to  extensive  melting, 
or  significant  expansion  of  the  polar  cap,  respectively. 

Faegre  (19  72)  obtained  for  a  certain  given  set  of  values 
of  the  parameters  five  distinct  solutions  of  his  variant  of  the 
Sellers  time-independent  model.   Two  of  these  were  highly 
asymmetric,  and  disappeared  when  c(()),u)  in  (1)  was  taken 
as  constant;  hence  Faegre  considered  these  solutions  to  be 
spurious  and  unphysical.   It  was  the  desire  to  eliminate  a 
priori  such  solutions  that  suggested  the  choice  of  symmetric 
coefficients.   Faegre 's  formulation  of  the  albedo  is  slightly 
different  from  that  of  Sellers,  mainly  in  that  the  minus 
subscript  in  (2a)  (i.e.,  the  cutoff  of  a((|),u)  at  u  )  was 
missing. 

The  three  symmetric  solutions  of  Faegre   could  be 
described  as  corresponding  to  the  present  climate,  an 
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ice-age  climate  (about  15°  C  colder  on  the  average  than  the 
previous  one),  and  a  completely  ice~covered  earth  (at  an 
average  temperature  of  about  175  K) .    This  last  climate 
was  also  the  one  obtained  by  Sellers  when  decreasing  the 
solar  radiation  by  more  than  4%. 

These  results  raised  the  question  of  the  transitivity 
of  the  earth's  climate,  as  formulated  by  Lorenz  (1968,  1970). 
In  Lorenz' s  terminology,  a  time-dependent  system  of  equations 
is  transitive  if  its  solutions  have  a  unique  equilibrium 
statistic,  that  is,  if  all  solutions,  independently  of 
initial  conditions,  have  the  same  infinite  time  average; 
otherwise  the  system  is  intransitive.   Lorenz  pointed  to  tne 
existence  of  certain  transitive  systems  which  possess  a 
property  called  by  him  almost  intransitivity ,  i.e.,  that 
of  having  at  least  some  solutions  whose  averages  over  long, 
but  finite,  time  intervals  are  different  —  these  solutions 
then  would  alternate  in  time  between  the  different  averages. 
He  raised  the  possibility  that  the  atmospheric  system  is 
almost  intransitive;  in  other  words,  that  the  known  climate 
changes  in  the  past  were  not  necessarily  caused  by  changes 
in  external  conditions  (like  solar  radiation),  but  rather 
were  an  effect  of  the  normal  evolution  of  the  system. 

Schneider  and  Gal-Chen  (197  3)  investigated  the  question 
of  transitivity  for  energy-balance  climate  models  by 
formulating  time-dependent  versions  of  the  Budyko  (B) , 
Sellers  (S  and  SV)  and  Faegre  (F)  models.   They  solved 
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numerically  the  initial-value  problem  governing  these  time 
dependent  models  for  a  large  range  of  initial  conditions. 
The  models  were  found  to  be  intransitive,  rather  than  almost 
intransitive:   every  solution  tended  as   t  ^  °°   to  one  of 
two  (or  more,  in  the  case  of  the  (F)  model)  equilibrium 
solutions;  viz.,  the  equilibrium  statistic  of  the  system 
governing  each  model  was  not  unique,  and  no  spontaneous 
transition  from  one  equilibrium  to  another  was  possible. 
The  two  equilibrium  solutions  obtained  for  all  models 
corresponded  to  the  present  climate  and  to  the  previously 
mentioned  completely  ice-covered  earth. 

The   equilibrium  solutions,  at  least  for  the  (S)  and 
(SV)  models,  proved  stable  under  rather  large  perturbations 
in  both  initial  conditions  and  parameters.   That  is,  solu- 
tions which  differed  in  their  initial  conditions  from  one 
of  the  limiting  "equilibria"  by  as  much  as  +  18  K  tended 
as  t  ->•  «>   to  the  "equilibrium"  near  which  they  started. 
Also  changes  of  +  1.5%  in  the  solar  radiation  led  to  limit- 
ing equilibria  which  differed  by  only  a  few  degrees  from 
the  unperturbed  ones.   However,  changes  of  more  than  -  18  K 
in  initial  conditions  or  -  1.5%  in  solar  radiation  led  from 
the  "present  climate"  to  the  "ice-covered  earth".   The  latter 
seemed  to  be  the  most  stable  climate  in  all  investigations 
mentioned  above . 

Schneider  and  Gal-Chen  did  not  obtain  a  limiting  steady- 
state  solution  which  would  correspond  to  a  true  ice  age  as 
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recorded  in  the  planet's  history.  Their  results  with  the 
time-dependent  version  of  the  Faegre  model  (F)  were  rather 
different  from  those  with  the  two  versions  of  the  Sellers 
model  (S  and  SV) ,  especially  with  regard  to  the  stability 
of  steady  states  under  perturbations. 

Contrary  to  the  results  of  Schneider  and  Gal-Chen, 
Dwyer  and  Petersen  (1973),  with  a  time "dependent  model 
essentially  identical  to  Schneider  and  Gal-Chen's  (S) , 
obtained  solely  one  type  of  limiting  steady  state,  that 
corresponding  to  the  "present  climate".   They  used  only  the 
data  climate   u  =  uC^))   as  initial  conditions,  but  varied 
the  solar  radiation  Q.   The  actual  values  of  the  limiting 
equilibrium  depended  of  course  on  the  values  of  Q  used, 
but  slight  changes  in  Q  yielded  only  a  difference  of  a  few 
degrees  between  the  average  temperature  of  the  equilibrium 
and  that  of  the  data  climate;  no  "deep  freeze"  equilibrium 
was  obtained. 

These  results  seemed  to  be  interesting  enough  in  order 
to  warrant  further  study  of  energy-balance  climate  models. 
As  indicated  in  the  previous  section  we  chose  to  investigate 
symmetric  solutions  of  the  (SV)  model  of  Schneider  and 
Gal-Chen,  and  of  some  variations  thereof,  including  the  (F) 
model.   We  hope  that  this  study  will  add  as  much  light  and 
as  little  heat  as  possible  to  the  climate  question  (Jackson, 
1962) . 
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4 .      Steady-State   Solutions 

We  turn  now  to  the  mathematical  theory  of  equation  (1) . 
Introducing  the  new  space  variable  x  =  2(})/tt,  we  obtain  the 
initial-and-boundary   value   problem    (IBVP) 


C(x)u^   =    (-)^       ■     /     ,.,,    |-  sin   ^-[k-  (x)+k^(x)g(u)  ]  u^ 
t  tt'       sin(iTx/2)     3x  2         1  2         ^  x 


-   ou    [1   -   m  tanh    (c^u   ) ] 
(4) 


+    MQ{x){l-b(x)+Cj^[Uj^+(u-C2z(x)-u^)_]  j-^    , 


0<x<l,       0<t, 


(5a)       u    (0,t)    =   u    (l,t)    =    0    ,       (5b)       u(x,0)    =    u(x)     , 


where    (4)    is    a  nonlinear  parabolic  partial   differential 
equation    (PDE) .      Here      g(u)    is    given   by    (2e) ,       y   =    1 
(its    significance  will    appear    later,    in   Section    6)    and 


c,    =   0.009,      c^   =   0.0065    deg   m"-"-,         0^=    1.9    x    10  deg      , 

3  -2 

(6)       c.    =    6.105x0. 75xexp(19. 6) xlO    dyn   deg  cm      ,    c^=    5350    deg, 

a      =    1.356xl0~"'-^cal    cm~^sec~"''deg~'^ ,    m  =    0.5,    u^=    283.16    deg, 


Mesh  data  at  10°  latitude  for  C(x),  kj^(x),  k2(x),  Q(x),  b(x), 
z(x)  are  given  in  Table  1.   The  units  of  the  constants  and 


-16- 


mesh  data  are  chosen  such  that  the  common  units  of  all  terms 

-2    -1 
in  (4)  are  cal  cm  sec 

.  The    first   step   of   the    investigation    is    to   find   steady- 
state   solutions    of    (4,    5a)    in    the    range   of   physical   interest, 
and    to  study    their   dependence    on   chemges    in   the   model.      We 
consider   therefore    the   steady-state   equation   obtained    from    (4) 
by   setting      u .     e    0 .      After   some    rearrangement  we   get   the 
following   two-point  boundary-value   problem    (BVP)     for   the 
system  of   ordinary   differential   equations    (ODE): 


(7a)  u      =    V 


V 


,.,2    F(x,u)  TT.  .X.  k;(x)+k'(x)g(u) 

k- (x) g'  (u)       2 

Vt ^ V       ,  0    <    x    <    1    , 

k(x,u) 

(8)  v(0)    =    v(l)    =    0    , 


where 


(7c)  k(x,u)    =    kj^(x)    +    k2(x)g(u)     , 

(7d)  F(x,u)    =    yQ(x)-(l    -   b(x)+    c^  [u^+     (U-C2Z  (x) -u^^)  _  ]  |^ 

4  6 

-   ou    [1-m   tanh    (c^u   ) ]    . 


We    use   shooting    (Isaacson    and  Keller,    1966,    Keller,    1968) 
as    the   numerical   procedure   to   solve    (7,    8) :    equations    (7a, b) 
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are  solved  with  initial  conditions 

(9a)    v(0)  =  0  ,  (9b)   u(0)  =  u^  , 

for  different  values  of  u^ ;   denote  by  u(x;u^),  v(x;u„) 
the  solution  of  the  initial-value  problem  (IVP)  (7,9)   in 
order  to  emphasize  its  dependence  on  the  parameter   u^.. 
An  iterative  scheme  is  then  used  to  obtain  those  values 
of   Uq  which  satisfy 

(10)  v(1;Uq)  =  0  . 

For  these  values  of   u^   the  solution   u(x;Uq),   v(x;u-) 

of  the  IVP  (7,9)  is  also  a  solution  of  the  BVP  (7,8). 

To  obtain  numerical  solutions  of  prescribed  accuracy  to 

the  BVP  (7,8)  one  has  therefore  (Keller,  1968)  to  achieve 

the  desired  accuracy  both  in 

(i)    solving  the  IVP  (7,9),  and  in 

(ii)   solving  iteratively  the  nonlinear  (finite)  equation  (10) 

In  solving  numerically  the  IVP  (7,9),  we  encounter  the 
difficulty  that  (7b)  is  singular  at  the  origin,  since 
cot(iTx/2)  ^  0°  as  X  ^  0 .   This  singularity  arises  from  the 
mathematical  form  of  the  diffusion  term  D  of  (1)  in  spherical 
coordinates.   It  is  easily  shown,  though,  that  (analytical) 
solutions  of  the  IVP  exist  and  are  unique  at  least  "in  the 
small",  i.e.,  for   | Uq-Uq |  <    c    ,      0  <  x  <  6,  arbitrary   U- 
and  some   6,  e  >  0.   The  numerical  difficulty  of  the  initial 
point  being  singular  can  be  circumvented  by  using  a  variable- 
step  method. 
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It  also  turns  out  that  in  order  to  achieve  prescribed 
over-all  accuracy  of  the  numerical  solution   of  (7,9)  with 
a  minimum  of  computational  effort  it  is  convenient  to  use  a 
variable-step  variable-order   multistep  scheme.   The  ODE 
solver  we  computed  with  was  developed  at  the  Lawrence  Liver- 
more  Laboratory  and  is  documented  by  Hindmarsh  (1972,1974). 
Both  the  Gear  and  Adams  methods  included  in  the  package  were 

used  and  gave  results  identical  to  within  the  prescribed 

-7 
relative  accuracy,  which  was  10   ;  however,  the  Adams  method 

was  faster  and  was  therefore  used  in  most  integrations. 

The  number  of  steps  needed  per  solution  (per  value  of  u„) 
for  this  accuracy  was  of  the  order  of  100.    It  is  of  interest 
to  point  out  that  near  the  singularity  at   x  =  0,   and  near 
the  point  where  the  jump  discontinuity  in  the  derivative  of 
F(x, u(x; Uq) )  occurred,  the  order  of  accuracy  chosen  by  the 
scheme  was  one  and  the  step  size  was  very  small,  so  that  a 
proportionately  larger  amount  of  computation  was  done  in  the 
neighborhood  of  these  two  points. 

Thus  every  solution  of  (7,9)  for  given    u   puts  at 

our  disposal  a  value  of  the  function   v(1;Uq),  accurate  to 

-4  * 
10   .    To  find  the  zeros  of  v(1;Uq) ,  we  used  the  method  of 

false  position  (Isaacson  and  Keller,  1966).   The  criterion 

for   u„   to  be  a  root  of  (10)  was  |v(l;u„) |  <  10 


In  the  units  chosen,  typical  values  of  u  are  0(10  )  and 
of  V  are  O(IO^) . 
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The  range  of  physical  interest  in  which  we  searched  for 
solutions  of  (10)  was  lOOK  £  u^  _<  300K  (see  also   Faegre, 
1972).   Within  this  range,  three  solutions  of  the  BVP  (7,8) 
were  found  which  are  denoted  by   u, (x) ,  u„(x),  u^(x). 
They  correspond  to 

u,(0)  =  247. 74K,  U2(0)  =  223. 97K  ,  ^^^^^^    "  168.94  K. 

The  curve  u  (1)  vs.  u(0),  the  zeros  of  which  correspond 
to  the  solutions   Uo,U2,u,  ,  is  given  in  Figure  2.    The 
individual  solutions   u, (x) ,  U2(x),  u^(x)   are  plotted  in 
Figure  3. 

It  is  customary  and  useful  to  characterize  a  climate 
of  the  model,  i.e.,  one  of  the  solutions  above,  by  its 
average  temperature,  rather  than  by  the  temperature  at  the 
pole.   Therefore  we  introduce   for  functions   cf)   the 
averaging  operator   £   by 

2  ,2 


£4.  =  J  (j)(x)  sin  {—)    dx/  J 
n  /   n 


sm  (~y)    dx  . 


0  ^    0 

In  particular,  it  is  known  (Isaacson  and  Keller,  1966)  that 

for  functions  (j)  symmetric  about  the  equator,   (|)(2-x)  =  (p  {x)  , 

which  also  satisfy   4)  (0)  =  0,  and  hence  can  be  extended  so 

as  to  have  period  2,  the  trapezoidal-quadrature  approximation 

of   S   is  very  accurate.   Denote  this  approximation  by  &      , 

where   A  =  1/9  (corresponding  to  a  10°  mesh) ;  then  we  have 


S,u,  =  287.76  K,   £.u„  =  267.44  K,   S.u-.  =  175.43  K  . 
A  1  A  2  'A3 
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Clearly  u.  corresponds  to  an  interglacial,  u„  to  a  glacial, 
and  u^  to  the  ice- covered  earth. 

Note  that  for  the  data  climate   u  =  u(x) 


u(0)  =  247.36  K,    £^u  =  287.20  K  , 

which  is  indeed  very  close  to  the  interglacial  solution  of 
the  model,   u-,  (see  also  Figure  3).   Though  present  observa- 
tions were  used  to  obtain  constants  and  empirical  functions 
in  (4) ,  no  explicit  term  was  added  to  ensure  that  (4)  hold 
with   u   =  0;  also  the  diffusion  term,  D,  is  of  the  same 
order  of  magnitude  as  the  radiation  terms,  R.  and  R-  ,  so 
that  this  is  not  the  result  of  a  simple  radiation  balance. 

The  same  is  true  of  the  work  of  all  the  previous  authors 

* 

discussed   in   Section    3    ;    these    authors,    however,    assumed 

at   least   implicitly      that   the   data   climate   should   be    a 
steady-state   solution   of    the   model,    or   approximately   so. 
Since    this    result   is   not    actually  built   into   the   model,    as 
far    as   we   can   see,    we    think    it   is    rather   remarkable    that 
this    class   of   models   yields    a   steady-state    solution,    u, (x), 
very    close    to   the    data    climate,    u(x).    Henceforth  we    shall 
refer   to   u,    also   as    the    "present   climate"    of    the   model. 


* 

In  fact,  Schneider  and  Gal-Chen  (1973)  did  use  an  extra 

"fudge  coefficient",  0.97  £  Cq(x)  <    1.03,  in  R^  in  order 
to  achieve  better  agreement  between  the  interglacial,  or 
"present"  climate  of  their  models  and  the  data  climate; 
they  obtained   £u,  =  287.06  K  vs.  £u  =  287.30  K  for  the 
( S )  mode 1 . 
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The  mesh  data,  from  which  constants  and  empirical 
functions  for  the  model  were  computed,  are  known  to  be  in 
error  by  possibly  as  much  as  100%  (Schneider  and  Gal-Chen, 
1973);  also  some  of  the  parameterizations  used  in  formulating 
the  model  are  questionable.   Therefore  we  made  a  number  of 
computations  in  order  to  obtain  information  on  the  dependence 
of  the  physically  significant  solutions,  u, ,Up,u^,  on  varia- 
tions in  the  model.   The  salient  features  of  these  computa- 
tions, summarized  in  Table  2,  are: 

(1)  The  dependence  of  all  these  solutions  on  the 
coefficients  seems  to  be  very  smooth. 

(2)  The  bounds  on  the  albedo,  0.25  £  a(x,u)  <_   0.85,  are 
essential  for  the  existence  of  three  steady-state  solutions 
in  the  physical  range,  100  K  £  u(0)  <_   300  K:  u^  disappears 
when  the  bound  a  <_  0.85  is  not  enforced,  and  u,  disappears 
when  the  bound  0.25  £  a   is  not  enforced. 

(3)  The  terrestrial  radiation  term  Rq  ,  given  by  (Ic), 
and  its  nonlinearity  are  essential  for  the  existence  of 
physically  significant  solutions:   if  the  term  is  set  equal 
to  zero  all  solutions,  u, ,U2/U   ,  disappear;  when  it  is 
replaced  by  its  average, 

R-  =  £  c(x,u(x))au  (x)  =  const.  , 

U2  only  obtains;  when   Rq  is  replaced  by  a  linear  approximation, 

Rq  =  au-^{u+4  (u-u)  }£.c(x,u(x)  )  ,    u  =  £^u(x)  , 
then  u,  only  obtains. 

-22- 


(4)  Our  results  with  the  model  using  the  Faegre  formula- 
tion of  the  albedo  are  very  similar  to  those  using  the  Sellers 
formulation;  it  is  hard  to  explain  the  disagreement  in  this 
respect  between  our  results  and  those  of  Schneider  and 
Gal-Chen  (1973) ,  who,  as  mentioned  in  Section  3,  obtained  very 
different  results  when  using  the  two  albedo  formulations. 

(5)  The  values  of  the  empirical  functions  in  the  model 
when  fitting  mesh  data  by  (i)  Bernstein  polynomial  approxima- 
tion and  (ii)  cubic  spline  interpolation  are  pointwise  rather 
different  for  some  of  the  functions  (see  Figure  1) .   The 
solutions  of  (7,8)  when  using  these  different  fits  are,  however, 
very  close.   This  also  supports  the  assertion  in  paragraph  (1) 
above,  and  shows  that  the  uncertainty  in  the  mesh  data  does 

not  affect  in  an  important  way  the  conclusions  of  the  investi- 
gation. 

We  explored  another  variant  of  the  model,  in  which  the 
eddy  diffusivity  k(x,u)   given  in  (7c),  corresponding  to  the 
(SV)  model  of  Schneider  and  Gal-Chen,  was  replaced  by  k(x,u), 
where   u   is  the  data  climate;  hence  also  g'(u)v  in  (7b)  is 
replaced  by  g'(u)u'(x).   This  variant,  in  accordance  with  the 
terminology  of  Gal-Chen  and  Schneider  (1975)  ,  would  correspond 
to  an  (SVC)  model.   We  denote  this  modification  of  (7)  by  (7*); 
the  solutions  of  (7 ',8)  are 

u^(0)  =  247.55  K,  U2(0)  =  227.76k,  U3(0)  =  169.44  K, 
£  u   =  287.70  K,   £  U2  =  268.60  K,   £^U3  =  175.44  K. 
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These  are  indeed  very  close  to  the  solutions  of  (7,8), 
especially  for  u,  and  u^  (see  also  Figure  3)  . 

Gal-Chen  and  Schneider  (19  75)  investigated  the  effect  of 
variations  in  the  solar  radiation  yQ  on  the  equator-to-pole 
temperature  gradient;  they  argue  that  this  dependence  should  be 
monotone,  in  fact  monotone  increasing.   In  the  light  of  this 
argument  and  of  their  results,  model  (SV)  is  superior  to  (S), 
as  already  mentioned  in  Section  2;  moreover,  model  (SVC)  is 
superior  to  (SV) .  Since   (SVC)  is  also  more  convenient  to  use 
when  investigating  the  stability  of  the  steady-state  solutions, 
we  shall  work  with  it  in  the  sequel;  the  corresponding  form 
of  (4)  we  denote  by  (4')^  the  corresponding  form  of  (1)  by  (1'), 
and  the  corresponding  form  of  (7)  by  (7'). 

For  this  model,  additional  computations  were  made  outside 
the  range  of  physical  interest.   One  more  steady-state  solution 
was  found;  we  denote  it  by  u. (x),  and  it  satisfies 

u^(0)  =  -  185.99  K,    £^u^  =  -  175.40  K  . 

Most  probably  there  are  no  other  solutions  of  the  BVP  (7', 8) 
at  all.   Indeed,  the  solutions  of  the  IVP  (7', 9)  become 
unbounded  as   u„   moves  towards  the  ends  of  the  range  explored, 
viz.,  -1330  K  £  Uq  £  300  K;  i.e.,  u(x;Uq)  ^  +°°  as  Uq  approaches 
the  ends  of  the  interval  above  (see  Figure  2). 

We  would  like  also  to  point  to  the  fact  that  some  of 
the  modifications  of  the  model  which  have  three  physical 
steady  states  yielded  numerical  values  of  the  latter  consider- 
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ably  higher  than  u,  ,u_,u^  (see  Table  2);  however,  it  was  the 
qualitative  feature  of  the  number  of  solutions  and  their 
approximate  position  with  respect  to  each  other  that  was  of 
interest  in  our  investigation:   the  averaae  temperature  of 
a  given  solution  of  any  fixed  model  could  be  changed  to 
practically  coincide  with  that  of  the  solution  u.  of  (4)  to 
which  it  corresponds  by  adjusting  the  values  of  the  coeffi- 
cients. 
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5.      Stability   of   the   Steady-State   Solutions 

In    the   previous    section  we  have   shown    that  equation    (4') 
with    the   boundary    conditions    (5a)    has    three    steady-state 
solutions    of   physical   interest,    u   =   u.(x)  »    1   f.  j    £   3.       In 
this    section  we    shall   study    the   stability   of   these    steady 
states. 

Let  us  concentrate  on  any  one  of  the  steady  states  above, 
u  =  u.(x),  where  j  is  fixed.  Stability  of  u.  means  that  small 
perturbations    of   u.    die    out  with    time.      More      precisely,    u. 

o 

is  stable   if,  when  taking  any  nearby  state   U,  i.e.,  one 
which  differs  little  from  u.  , 

(11)  U  =  u. (x)  +  ev(x)  , 

o 

say,  where  v  is  arbitrary  and   e  >  0   is  not  too  l^rge,  then 
the  solution  U(x,t;e)  of  (4' -5)   with  initial  condition 

o 

U(x,0;£)  =  U(x;e)   tends  to   u.  itself  as   t  >  «=  . 
Let  (4')  be  written  symbolically  as 

(12)  u^  =  N(u)  , 

where  we  divided  through  equation  (4')  by  C (x) ,  0  <  C  £  C (x) 
_<  C„  <  »,  and  N(u)  is  the  corresponding  right-hand  side.  Take  a 
perturbed   u,   u  =  U(x,t,e:),   which  is  assumed  to  satisfy 
the  boundary  conditions  (5a)  and  the  initial  condition  (11) . 
If  such  a  solution  of  (12)  exists  for  all  e  sufficiently  small, 
0  ^  ^  £  ^Q  '  then  the  equations 
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U  (x,t;e)-U. (x,t;0)      N(U(x,t;£))-N(U(x,t;0)) 
(12')   ^ S =    


0  <  e  ^  Eq  ' 

also  hold,  with  U(x,t;0)  =  u.(x).   Letting    e  ^  0   we 
obtain  the  linearized  equation 


(13a)  V  =  -  L  .V  , 

where  we  define   v(x,t)  s  -r-^  U(x,t;£)|  _„  and 


(13b)   L   =  -  |^N(u(x,t;e))|^^Q  =    -   |^N(u)|^^^   ,  l<j<3, 


At  this  point,  for  the  sake  of  simplicity,  we  shall 
drop  the  subscript  j,  and   u  will  thus  stand  for  some   u.  , 
L   for  the  corresponding  L . . 

Equation  (13a)  is  linear  in  v  and  it  has  a  unique  solu- 
tion  V   satisfying  the  boundary  conditions 

(13c)  V  (0,t)  =  V  (l,t)  =  0  , 

and  arbitrary  initial  conditions 

(13d)  v(x,0)  =  v(x)  . 

The  solution  may  be  found  by  the  usual  method  of  separation 
of  variables,  or  expansion  in  eigenfunctions  (normal  modes). 
Consider 


v(x,t)  =  e    w(x)  ; 
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this  V  will  be  a  solution  of  (13)  with  v(x,0)  =  w(x)  iff 
(if  and  only  if)  X  is  an  eigenvalue  and  w  is  an  eigenf unction 
of  L,  i.e.,  iff  the  pair  (A,w)  satisfies  the  homogeneous 
problem 

(14a)  Lw  =  Aw 

(14b)  w  (0)  =  w  (1)  =  0  . 

We  shall  show  that  the  operator  L  has  sufficiently  many 
eigenfunctions ,  and  that  therefore  any  solution  v  of  (13) 
can  be  written  as  a  series 


(15)  v(x,t)  =  I      a.  e~^    "^   w^^^  (x) 

k=0   ^ 

(]c)  (]r) 

where  (\         ,w^  ')  are  all  the  solutions  of  (14).   It  is 
clear  from  the  derivation  (12')  of  (13)  and  from  (15)  that 

for  the  steady  state   u  to  be  stable  it  is  necessary  that 

(k) 
every  eigenvalue  X  of  L  have  positive  real  part: 

(16)  Re  A^'^^  >  0  . 


This  condition  defines  the  so-called  linear  stability  of  u. 

It  has  been  shown  rigorously  in  a  few  cases  and  it  is  believed 

in  many  others  that  (16)  is  also  sufficient  for  the  stability 

of  u  as  we  defined  it  at  the  beginning  of  this  section,  i.e., 

that  linear  stability  implies  nonlinear  stability.   In  the 

next  section  we  shall  give  an  argument  which  will  make  it  at 

least  plausible  that  this  is  the  case  also  for  the  problem 

at  hand. 
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We  turn  now  to  the  analysis  of  the  eigenvalues  of  the 
linear  second~order  ordinary  differential  operator  L,  which 
we  can  write  as 


(17a) 


Lw  = 


^    [p(x)wj^  +  q(x) 


w 


Here  p,q,r   are  determined  by  (13b)  and  by  (4'): 


(17b) 


and 
(17f) 


TTX, 


r(x)  =  C(x)  sin  (~2~)  ' 


(17c)         p(x)  =  (|-)^  sin  (^)-k(x,u)  , 


(17d) 
with 

(17e)    (c 


q(x)  -  {Q(x}{c^)^,    -    d(x,u)}/C(x)  , 


I'c' 


f    c,   if  0.25  <_   a(x,u)  <_   0.85 
-    \  and  u  (x) -CpZ  (x) -u  <_  0    , 

*■  0    otherwise. 


d(x,u)  =  T—  c(x,u)au 

oU 


=   au   <4[l-m   tanh(CoU    )]-    Smc^u    [1-tanh    (c^u    )]>, 


Clearly   L    thus    defined   is    formally    self-adjoint    (Courant 
and   Hilbert,    1953,      Birkhoff    and   Rota,    1969)       under    the 
prescribed  boundary   conditions   with    respect  to   the    inner 
product 

(18) 


(w-,,W2)    =         r(x)w    (x)w2(x)    dx 


Moreover,    r,p      are    nonnegative    and    twice    continuously   differ- 
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entiable  on  the  interval  I  =  [0,1],  and  q  is  piecewise  contin- 
uous on  I.  However,  because  of  the  singularity  at  x  =  0  due  to 
the  fact  that  r(0)  =  p(0)  =  0, the  usual  Sturm- Li ouvi lie  theory 
of  self-adjoint  operators  (Courant  and  Hilbert,  1953,  Birkhoff 
and  Rota,  1969)  does  not  apply  to  L;  this  difficulty  though 
can  be  overcome  and  the  theory  can  be  extended.   The  crucial 
element  in  this  extension  is  the  observation  that,  because 
of  the  boundedness  of  q,  L  is  bounded  from  below  in  the  sense 
that,  for  any  w  satisfying  (14b)  for  which  (w,w)  <  °° ,    the 
inequality 

(19a)  <w,w>  >^  K(w,w)  , 

holds  with  some  fixed  constant  K,  K  >   min  q(x)  >-°°,  independent 

0£xj<l 
of  w;  here  <T'r,w>  is  the  Dirichlet  integral  corresponding  to  L, 

1 

f  2  2 

(19b)       <w,w>  =  (Lw,w)  =    (pw   +  rqw  )  dx  . 

0 

This  result,  together  with  an  analysis  of  the  nature  of 
the  singularity  at  x  =  0,  are  sufficient  to  show  that  indeed 
L  has  a  complete  system  of  eigenf unctions ,  orthonormal  with 
respect  to  the  inner  product  (18)  (Birkhoff  and  Rota,  1969) , 
and  that  the  eigenvalues  of  L  are  real  and  can  be  arranged 
in  an  ascending  sequence 

c=o  ^  >(1)  ^  i  (2)  ^  ,  (3)   . 

with  X^^^    ->-  <»   as   k  ->  «>,  and  K  =  X^"""^   in  (19).   Hence 
any  solution   v  of  (13)  can  be  written  as  a  series  (15) . 
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Therefore  the  stability  question  for  the  steady  state  u 
reduces  to  determining  whether 

(16')  A^^^  >  0  , 

in  which  case   u   is  stable,  or  whether  the  opposite  holds 
in  which  case   u   is  unstable. 
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5a.      Stability   Criterion 


In    this    subsection  we   shall   show   that    it   is   possible    to 
determine    the   sign   of    A  without   actually   computing    A         . 

We    start  with    the   variational   characterization   of   the 
lowest  eigenvalue    (Courant   and  Hilbert,    1953),    A         , 


(2  0)  A^^    =        min        ']^'^^   =        i"in         <v,v>    , 

<v,v><«>    ^^'^''  (v,v)=l 

(v,v)^0 

i.e.,   A     is  the  minimum  of  the  Rayleigh  quotient  R(v) 

corresponding  to  L, 


(21)  R(v)  =    <v,v>/(v,v)  , 

and  the  minimum  is  assumed  for  v  =  w    .   Indeed,  (14a)  is 
the  Euler  equation  for  (20) ,  and  v  (1)  =  0  is  the  "natural" 
boiandary  condition  for  (20)   in  the  sense  of  the  calculus  of 
variations  (Courant  and  Hilbert,  1953);  also   v  (0)  =  0  is, 
according  to  the  theory  of  singular  operators  with  the 
properties  of  L,  the  only  boundary  condition  at  x  =  0  which 
ensures  that,  for  solutions  v  of  (14),  <v,v>  as  well  as  (v,v) 
is  finite.   In  particular,  we  conclude  that  the  minimizing 
function,  v  =  w    ,  is  at  least  as  smooth  as  the  coefficients 
of  L,  viz.,  it  must  have  a  piecewise  continuous  second 
derivative. 

With  these  preliminaries  we  are  able  to  prove  the 
following  known 
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Lenma.    The  first  eigenfunction  of   L,   w    ,   is  strictly 


positive , 


w^-'-^x)  >  0  ,  0  <  X  <  1  . 


Proof:    From  the  definition  (21)  of    R(v)   it  is  clear  that 


R(  iw^-"-^  I  )  =  R(w^-'-^ 


/ 


where   |y|   denotes  the  absolute  value  of  y.   If  w     had  a 
zero  at  some  interior  point   >-'q  /   0  <  x.  <  1,  and  w    {x    )f^0, 
then  the  first  derivative  of  |w    |   would  have  a  jump  at 
X  =  Xq  ,  which  contradicts  the  smoothness  of  |w    |  as  a 
solution  of  the  variational  problem.   If,  on  the  other  hand, 
w^"""^  (Xq)  -  0,   0  <  Xq  <  1,   or   w^-"-^  (0)  =  0,   or  w^"""^  (1)  =  0, 
then,  by  the  uniqueness  theory  of  linear  ODE  (Birkhoff  and 
Rota,  1969),  it  would  follow  that  w     e  0;  this  contradicts 
the  fact  that  w     is  a  nontrivial  solution  of  (14) , 
completing  our  proof. 

Now  we  are  ready  to  state  our  stability  criterion,  which 
is  part  of  the  conclusion  of  the 

Theorem.   Let   L,   A=A    ,   w=w     be  as  above.   Suppose 
also  that  there  exists  a  function   v,   as  smooth  as  w, 
nonnegative,   v  >  0   on  I ,  and  satisfying 


(22)        Lv>0,    v(0)=0,    v(0)  =  l. 

Then   v  (1)  >^  0  implies  A  >  0.   Moreover, 
(i)    A  >  0  if  either   Lv  ^  0   or   v  (1)  >  0,   and 
(ii)   if  Lv  E  0  holds,  then  A>0,  A=0,  or  A<0,  according  to 
whether   v  (1)  >  0,   v  (1)  =  0 ,  or   v  (1)  <  0. 

X  '      X  X 
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Proof;         The    required   results    are   easily   read   off    from   the 
following   sequence    of   eq       .ities    obtained    from  the   definitions 
and  by    integration   by   parts: 


A        rvw  = 
0 


A(w,v)    -    (Lw,v)    =    I    -  (pw    )    V   +    rqwv 


1         1 


-pw   V 

'^     X 


=      pwv 


+    I    pw    V      +    rqwv 

2\,      2\ 


0  0 

1  1 


+    I    -  (pVj^)  jjW   +    rqvw 


0         0 


=       (pwv    ) (1)    +     (Lv,w)     , 


where  we    used        w    (0)    =  w    (1)    =    0      and      v    (0)    =    0    . 

X  X  X 

Note.    This  is  essentially  a  comparison   theorem,  of  the 
type  familiar  from  Sturm- Li ouvi lie  theory  (Birkhoff  and  Rota, 
1969)  . 

The  result  under  (ii)  above  is  the  stability  criterion 
which  we  shall  use.  For  completeness  we  give  here  also  the 
following  slight  generalization  as  a 


Corollary.   Let   L,  X,  w  be  as  in  the  Theorem.   Suppose  v 

satisfies  the  hypotheses  of  the  theorem,  except  for  that  of 

being  nonnegative  on  I.  Instead,  let   x.  ,  0  <  x^  _<  1 ,  be 
the  first  zero  of  v, 

v(x)  >  0  on  0    <_  X   <   Xq. 
Then   X  <  0. 
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Proof:    The  result  follows  from  the  same  integration  by 
parts  carried  out  in  the  proof  of  the  theorem,  except  that 
now  the  upper  limit  of  integration  has  to  be  x^^  rather  than 
1,  and  we  use  v(x„)  -    0  rather  than  w  (1)  =  0  in  passing 
from  the  second  to  the  third  line.   Moreover,  since  v(x)  >  0 
for  X  <  x„  ,  and  v  is  continuously  dif ferentiable ,  v  (Xq)£  0. 
Furthermore,  it  cannot  be  that  both  Lv  =  0  and  v  (x.)  =  0, 
since  then  we  would  have,  by  uniqueness,   v  =  0,  which 
contradicts  v(0)  =  1.   This  completes  the  proof. 
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5b.   Stability  Results 


In  this  subsection  we  shall  apply  the  stability  criterion 
obtained  in  Subsection  5a  to  the  steady -state  solutions  u.  , 
1  <  j  <  3. 

For  this  purpose  we  construct  functions  v.  by  solving 
(22)  with  the  equality  sign,  and  with  L  =  L.  given  by  {13b) 

and  (17).   The  results,  obtained  with  a  relative  numerical 

-7 
accuracy  of  10   ,  are  that 


(a)  V.  ,  j  =  1,2,3,  is  positive,    v.  ^  v.  >  0,    V.  >^  0.5  ; 

(b)  ^  v.(l)  >  0   for   j  -  1,3,   whereas   ^  v.(l)  <  0  for  j=2 
dx   3  J  '    '  dx   3  -■ 

The  actual  computed  values  are 

^  v^(l)  =  2.39970,   1^  V2(l)  =  -3.24312,  |3^V3  (1) -4.  31025 . 

These   values    of   the   derivatives    at   x  =    1,    together  with   the 
values    of  v.,    are   sufficiently   bounded   away    from   zero   in   order 
to  conclude    that    the    conditions   of   the    theorem  are   satisfied 
beyond   the   doubt   of  numerical   uncertainty. 

Hence    the   solutions    representing   the   present   climate    and 
the    ice-covered  earth    are   stable,    whereas    the   so-called   ice 
age    of   the   model    is    unstable.         This    result   agrees   with      and 
throws    additional    light   on    the    results    of  previous    authors,    in 
particular   the    time -dependent   integrations    of   Schneider    and 
Gal-Chen    (1973) . 
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Actually  the  stability  computations  were  carried  out 
also  for  the  solutions  of 


(4")  N'  (u.)  =  0  ,  1  1  J  1  3, 

with 
(23)     N'(u)  =  (-)   ^_^_-^_  sin  (^)-u^ 

+  Q(x){l  -  Bq+  C^uj^  -  CgOU^ 

where    the    siibscript   \     }        is    defined    in    (2c)    and 


Kq    =    S^k(x,u(x))    =    2.2x10    ^,    Bq=    g^b(x)    =    2.85881, 
Cg    =    £.c(x,u(x))    -    0.61    ; 


the  corresponding  linearization  of  N'  is 

(13-)         l'.    =    -   N     (u.)    =    (^)2       ■    y     ,.,    f-  sin(^)    \- 
^         '  j  u      j  tt'       sin(7Tx/2)     9x  *■    2'     8x 

+  Q(x)  (Cj_)^„  -  4cgau   , 
where 


(ci),..  = 


c^       if   0.25  <  1-B-+C,u.  <  0.85, 
1  —     0   1  J  — 

0    otherwise 


This  is,  according  to  the  experiments  carried  out  in  Section  4 

(see  Table  2) ,  the  simplest  form  of  (4)  which  still  yields 

approximately  the  same  steady -state  solutions  in  the  physical 

range  of  interest.   The  results  are  the  same,  i.e.,  u   and  u^ 
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are  stable,  and  u^  is  unstable.  It  seems  therefore  quite 
plausible  to  conclude  also  that  for  all  raodels,  lying  in  some 
sense  between  (4)  and  (4"),  the  steady  states  corresponding 
roughly  to  u, /U^jU^  have  the  same  stability  properties  as  the 
latter.  Combining  this  remark  with  the  one  made  at  the  end  of 
Section  4,  it  seems  that  we  have  a  result  about  the  stability 
of  the  steady  states  for  a  certain  type  of  energy-balance  model, 
rather  than  just  for  one  specific  model  of  this  type.  It  seems 
desirable  to  define  precisely  the  type  of  model  having  these 
properties,  and  we  intend  to  try  to  do  so  in  further  work. 

Having  thus  determined  the  stability  properties  of  the 
steady  states  of  (4'),  we  want  to  obtain  some  additional  infor- 
mation by  actually  computing  the  lowest  eigenvalues  X  .       ,  and 
corresponding  eigenfunctions  w.  ,  of  L .  ,  1  £  j  £  3.  The  eigen- 
values are 
AJ-'-^=  3.95390x10"^,  X^^^  =    -5  .  87662xlo"^  ,  X^^^^  =    5  .  13587xlo"^  , 

and  the  eigenfunctions  w.   ,  1  £  j  £  3,  are  plotted  in  Figure  4. 
The  method  for  computing    (A.   ,w .   )  was  again  shooting, 
this  time  with  respect  to  the  parameter  A.  That  is,  equation  (14a) 
with  L  =  L.  ,  was  solved  for  w{x;A)  with  initial  conditions 

w(0;A)  =  1  ,    w  (0;A)  =  0  , 
and  with  different  values  of  the  'shooting  parameter"  A.  The  zeros 

of  the  function  v  (1;A)  were  then  computed  by  regula  falsi  to 

-4 
within  an  accuracy  of  10   .  The  over-all  relative  error  in 

-7 
computing  solutions  of  the  initial-value  problem  was  10   , 

as  before. 
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Notice    from    (15)    that 

T  =  1/1  a]^^ 

is    a   characteristic   decay   or   relaxation    time    for   the   solution 
U(x,t;e)    of    (12)    to   u   =   u.(x).      This    time    is   of   the   order 
of   10   years    for   all    j, 

T-    =    8.0   years    ,         Tj    =    5.4   years,,         x^   =    6.2    years. 

In   this    context   it   is    remarkable    that   Schneider   and   Gal-Chen 
(1973)    state    that,    in   one   of   their   integrations,    the   solution 

o 

of  (4)  with  initial  data   u(x)  =  u(x)   and   y  =  0.984,  after 
dropping  rapidly  by  about  12  K  in  average  temperature,  was 
nearly  constant  for  about  50  years  of  simulated  time,  and 
then  finally  dropped  to  a  steady  state  close  to  our  U2(x). 
From  our  results  it  becomes  clear  that  the  solution  mentioned 
above  hovered  for  a  time  of  the  order  of  t^    around  U2  /  but 

could  not  persist  there  indefinitely  because  of  the  instability 

* 

of  Up  ,  and  finally  attained  u^  ,  which  was  stable. 

It  also  follows  from  (13)  and  (14)  that  multiplying  C (x) 
by  a  constant   k  >  0   will  result  in  t . ,   1  1  J  1  3,  being 
multiplied  by  k.   a  similar  statement  holds  also  for  the 
nonlinear  equation  (4) ,  since  such  a  constant   k   just 
corresponds  to  a  different  scaling  of  the  time  t  (see  also 
Schneider  and  Gal-Chen,  1973).   Experiments  in  determining 


We  shall  see  in  the  next  section  that  t^    increases  with 
decreasing  y . 
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characteristic   response    times    for   solutions    of    (4')    were 
done   by   Dwyer   and   Petersen    (1973),    who   used   two  heat   capaci- 
ties     C(x)  ,      both    larger    than    that   used  by      Schneider   and 
Gal-Chen    (1973)       and  by    us.       It   seems,    however,    that   upon 
decreasing      p      in    (4)    to      y   =    0.98,    as    they   always    took 

o 

u(x)  =  u(x)  in  (5b),  the  average  temperature  of  the  solution 

u(x,t)   dropped  rapidly  at  first,  and  only  slowly  thereafter, 

as  indicated  by  Figure  2  in  their  article.   Apparently  this 

slow  decrease,  which  shows  that  their  solution   u(x,t)  of  (4') 

was  approaching  a  steady  state  close  to  our   u„(x),  was 

interpreted  by  them  as  proving  the  nonexistence  of  a  steady 

state   u-,(x)  ,  which  had  been  obtained  in  the  work  of  Budyko, 

Sellers  and  Faegre   when   decreasing   y   by  a  similar  amount. 

The   influence    of   changes    in      y      on   the   steady-state 

solutions      u.    ,    1  £   j    £   3,      of    (4')    will   be   investigated   in 

the    following   section.      At   this   point  we   only  want   to   mention 

that,    whereas    a  multiplicative    factor      k      in      C(x)       affects 

the   magnitude   of      X  .  and  hence    of      x .    ,      our   theorem  shows 

D  3 

that  it  does  not  affect  the  actual  stability  of   u.  ,  i.e., 
the  sign  of   A  . 
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6.      Perturbed   Steady-State   Solutions    and  Bifurcation 

.    It   is    clear   that   steady-state    solutions    of    (4')    could 
not   spontaneously   evolve    into  each   other.      More   precisely, 
the    solution   of    (4'-5a)    with   initial    condition   u(x,0)    =   u.(x), 
j    =    1,2,3,       is      u(x,t)     =   u.(x).      Moreover,    it   stands    to 

o 

reason  that,  for  any  physical  initial  condition,  u(x)  >^  100  K, 
we  would  have 


lim  u(x,t)  =  u^ (x)  , 
t- 


--XO  -I 


with  j  =  1  or  j  =  3,  since   u~ (x)  is  (linearly)  unstable, 
whereas   u-,(x),  u^(x)   are  (linearly)  stable  (and  u^(x)  _<-170K 
<  0)  . 

Thus,  to  explain  physical  ice  ages,  one  has  to  consider 
perturbations  in  the  parameters  appearing  in  equation  (4'). 
Such  perturbations  would  presumably  be  caused  by  physical 
mechanisms  not  included  in  the  model.   The  most  likely 
candidate  for  such  a  role  is   p,   which  up  to  now  was  taken 
to  be  unity.   Indeed,  many  ice-age  theories  rely  heavily  on 
a  change,  however  small,  in  the  amount  of  solar  radiation 
reaching  the  lower  layers  of  the  atmosphere  (SMIC,  19  71, 
Beckinsale,  1965) . 

Some  attribute  the  assumed  decreases  in  solar  radiation 
to  changes  in  the  parameters  of  the  motions  of  our  planet 
(Milankovitch,  1930) ,   others  to  airborne   volcanic  dust 
due  to  an  increase  in  volcanic  activity  (Fuchs  &  Patterson, 
1947),  and  so  on.   There  has  also  been  concern  about  a 

-41- 


possible  climatic  catastrophe   being  imminent  because  of  the 
increase  in  the  quantity  of  industrial  pollutants  in  the 
atmosphere  (Rasool  and  Schneider,  1971)  . 

To  investigate  the  effect  of  such  changes  in  the  model 
at  hand,  the  curve   u  (1)  vs.  u(0)  was  recomputed  for 
different  values  of   y ,   in  particular  with  a  view  to  obtain- 
ing  u,  (x;y),  Up(x;ij).   One  important  result  is  that  these 
two  solutions  coallesce  for 

y   =  0.98152 

c 

and  disappear  entirely  for   y  <  y  .   The  bifurcating  solution 

u  =    u    (x)    =    u^(x;y    )    =   u^(x;y    )    was    computed  with   over-all 
C  X  c  ^  c 

-9       I  •    I        -7 
relative  accuracy  10    and   |u  (1)|  <  5xio   .    The 

(u(0) ,u  (l))-curve  corresponding  to  y  =  y   is  given  in  Figure  5; 

notice  that  it  is  very  flat  near  u(0)  =  u  (0),  which  makes 

-^  c 

it  difficult  to  compute   u  (x)  accurately. 

Further  computations  were  carried  out  for   y   =  0.982, 
0.985,  1.01,  1.02,  1.03,  1.04  and  1.05.   The  results  of  these 
computations  are  summarized  in  Table  3  and  plotted  in  Figure  6. 
It  is  quite  interesting  that,  whereas  for   y  >_  1.0   the 
dependence  of  u.(0),   a  u.(x),  j  =  1,2,  on  y  is  almost  linear, 
this  dependence  is  definitely  quadratic  near   y  =  y   ;  the 
latter  is  in  good  agreement  with  the  general  theory  of 
bifurcation  for  nonlinear  parabolic  problems  (Hoppensteadt 
and  Gordon,  1975).   According  to  the  theory,  there  exists 
in  principle  the  possibility  that,  instead  of  disappearing 
at  y  =  y   ,  the  two  solution  branches  u,(x;y),  U2(x;y)  could 
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merge  into  one  periodic  solution   u, -(^/tj-y)  for  p  <  M  • 
This  possibility  was  not  borne  out,   however,  by  the  time- 
dependent  computations  of  Dwyer  and  Petersen  (1973),  and  of 
Gal-Chen  and  of  Schneider  (1973,  1975) ;  we  did  not  investi- 
gate it  further. 

Concerning  the  problem  of  the  pole- to- equator  tempera- 
ture gradient, 

Au.  =  u.  (1)  -  u.  (0)  ,        1  <  j  <  3, 
:    1       D  -   - 

as  discussed  by  Stone  (1973)  and  Gal-Chen  and  Schneider  (1975)  , 
the  curve  Au.  =  Au.(ij),  j  =  1,2,  is  particularly  interesting. 
We  notice  that 

(a)  the  values  of  Au,  lie  below  those  for  Au-  ; 

(b)  the  values  of  Au-,  are  monotonically  decreasing  with  y; 

(c)  the  values  of  Au-  have  a  maximum  for  y  somewhere  between 
y  =  0.9  85  and  y  =  0.99; 

(d)  the  dependence  of  both  Au,  and  Au-  on  y,  but  especially  that 
of  Au,  ,  is  very  steep  near  y  =  y  . 

Being  aware  of  the  limitations  of  the  model,  as  pointed  out 
in  Section  2,  we  do  not  want  to  make  extensive  comments 
concerning  these  results,  but  only  note  them  for  comparison 
with  the  results  of  other  models  and  for  further  study. 

We  also  studied  the  stability  of  the  solutions  u  =  u^(x) 
and  u  =  u.,(x;y  ).   Repeating  the  construction  indicated  at 
the  beginning  of  Subsection  5b  for  v  =  v^  and  v  =  v^  ,  where 
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L   =  -  |nN(u;y^)|^^^   ,       j  =  c,3. 


we  obtained  the  following  results; 


(a)  V.  ,  j  =  c,3,  is  positive,   v .  >^  V .  >  0,   V.  >_  0.75  ; 

(b)  1^  v^(l)  =  -  7.56xlo"^  ,   1^  v^(l)  =  4.23362. 

Clearly   u^(x;)j  )  is  still  stable,  in  fact  (d/dx)  v^  ( 1;  y  ) 
=  4.23   is  very  close  to  (d/dx) v  (1; 1. 0 )  =  4.31,  showing 
the  extremely  smooth  dependence  of  u-,(x;p)  on  \x. 

The  negativity  of  (d/dx)v  (1)  would  seem  to  point  to 
outright  instability  of   u  (x) ,   but  in  fact  its  small 
numerical  value  indicates  that,  within  the  accuracy  of  the 
computations,  it  is  actually  zero,  i.e.,   u  (x)  is  neutrally 
stable.   Indeed,  the  mathematical  theory  of  nonlinear 
problems  (Nirenberg,  19  74)  shows  that  for  a  bifurcation  to 
exist,  as  in  the  case  at  hand,  the  linearization 


c 
of  the  spatial,  time-independent,  operator  N  at  the  bifurca- 
tion point  (u  ,y  )   has  to  have  a  zero  eigenvalue.  This 
follows  from  the  infinite-dimensional  generalization  of  the 
one -dimensional  fact  that  a  function  can  be  inverted  in  a 
neighborhood  of  a  point  at  which  it  has  a  nonzero  derivative , i-e., 
that  it  has  a  unique  branch  near  such  a  point. 

We  also  computed  the  lowest  eigenvalues   A .    and 

corresponding  eigenf unction  w.    of   L . (y  )  for  j  =  c,3, 

D         1   c 
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by  the  shooting  method  described  in  Subsection  5b.   In  this 
computation,  the  results  on  v.  mentioned  before  were  valuable 
in  making  a  first  guess  for  the  shooting  parameter  A.  The 
computations  yielded 

A^-"-^  =  -  1.287xl0~"'-''"  ,    A^""-^  =  5.07265x10"^  . 
c  3 

(1)  -9 

Again   we    notice    that    A^       (y    )    =    5.07x10         is    very    close    to 

A^-"-^  (1.0)    =    5.13xl0~^,    whereas    |  A^"*"^  |    =    0.4xio~^    a|"'"Ni.O) 

-2   (1) 
=  0.2x10    A„   (1.0)  is  practically  zero,  as  it  has  to  be 

analytically . 

It  is  clear  by  the  continuity  of  A .   (p)  in  the  parameter 
y  that  for   y   <  y  £  1.05   we  have   a|  ^  >  0 ,   A^    >  0,  and 
A2    <  0,   so  that  the  interglacial  and  the  ice-covered  earth 
are  stable  for  the  entire  range  of  y  explored,  whereas  the 
glacial  is  unstable  for  the  same  range  of  y.   Furthermore  the 
ice-covered  earth  is  stable  also  for  smaller  y. 

There  is  one  further  point  of  view,  which,  while  illumi- 
nating the  significance  of  the  neutral  stability  of  u  (x), 
also  argues  for  our  linear  stability  analysis  being  sufficient 
for  concluding  on  nonlinear  stability  or  instability  of  the 
steady-state  solutions  of  (4")  corresponding  to  different 
values  of  y.   This  viewpoint  has  to  do  with  the  existence  of 
a  variational  principle  for  (4').   Indeed 

N(u;y)  =  0 

is    the    Euler-Lagrange    equation    for    the   extrema    of    the    functional 

1 

J(u;y)    =         \p^x   "*"    ^G(x,u)  y    dx    ; 
0 
here   p,    r    are    given  by    (17b,c)    and 

-45- 


G(x,u)  ==    F(x,Lo)  dio  , 

where  F(x,u)  is  defined  by  (7d). 

Clearly  the  stable  solutions  u, (x;l),  u^Cxpl)  corres- 
pond to  local  minima  of  J(u;l),  whereas  Up(x;l)  is  a  local 
maximum.   As  u,  (x;ij)  ,  which  is  a  minimiam  for  y  >  y   , 
coallesces  with  u^(x;y) ,  which  is  a  maximum  for  y  >  y   /at 
y  =  y   ,  a  saddle  point  u  =  u  (x)  obtains,  whereas  U2(x;y  ) 
is  still  a  minimum. 

This  variational  interpretation  makes  it  very  plausible 
that  solutions  u(x,t;y)  of  the  "flow" 

u^=N(u;y)  ,    y>y^, 

with  initial  conditions  near  u.(x;y),  that  is  at  a  finite 
but  small,  rather  than  infinitesimal,  distance  from  u.(x;y), 
j  =  1,2,3,  will  tend  as  t  ->  °°  towards  u.  if  u  .  is  a  minimum, 
i.e.,  j  =  1,3,  and  away  from  it  when  u .  is  a  maximum,  i.e., 
j  =  2.   Similarly,  for  y  =  y   ,  solutions  starting  near 
u_(x;y  )  will  still  converge  to  u-,  ,  but  solutions  starting 
near  u  (x) ,  though  they  may  hover  for  a  long  time  near  u   , 
since  t,,t„  ->  «=  as  y  ->•  y   ,  will  eventually  move  away  from  it, 
along  a  negative  slope  of  the  saddle,  and  finally  tend  towards 
the  absolute  minimum  u-,.   This  seems  a  rather  satisfactory, 
although  heuristic,  explanation  of  the  results  of  the  time- 
dependent  computations  of  Dwyer  and  Petersen  (19  73)  and  of 
Schneider  and  Gal-Chen  (19  73) ,  which  we  mentioned  alreay  in 
Section  5. 
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7.   Concluding  Remarks 

We  studied  the  zonally- and- vertically  averaged  energy- 
balance  climate  model  governed  by  equations  (1-3) ;  these 
equations  are  based  on  simple  parameterizations  of  albedo, 
greenhouse  effect   and  eddy  diffusion  of  heat  in  terms  of 
yearly  averaged  sea-level  temperature,  which  is  the  only 
dependent  variable  of  the  model. 

Three  positive  steady-state  solutions  of  the  model, 
symmetric  with  respect  to  the  equator,  were  found  by 
accurate  numerical  computations,  and  apparently  no  more 
such  solutions  exist.   These  steady  states  can  be  identified 
with  an  interglacial  climate,  approximating  very  well  the 
one  prevailing  presently  on  earth,  a  glacial  climate,  and 
a  climate  during  which  the  earth  would  be  completely  ice 
covered.   The  climates  obtained  were  only  slightly  changed 
when  making  small  changes  in  the  numerical  values  of  the 
coefficients  and  when  making  certain  changes  in  the  functional 
form  of  the  model's  equations.   However,  the  bounds  on  the 
values  the  albedo  can  take  were  essential  in  order  to  obtain 
these  three  climates;  also  linearizing  the  outgoing  planetary 
radiation  resulted  in  a  reduction  of  the  number  of  solutions. 

We  then  determined  the  stability  of  the  time-dependent 

solutions  of  (1')  under  small  perturbations  about  the  model's 

steady  states.   This  stability  was  shown  to  depend  on  certain 

properties  of  a  comparison  function,  which  was  constructed 

numerically. 
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We  foiind  that  the  interglacial  and  the  "deep  freeze" 
climate  are  stable,  and  that  the  glacial  climate  is  unstable. 
This  means  that  the  first  two  can  obtain,  at  least  approximate- 
ly, as  steady  states  in  a  physical  system  governed  by  equa- 
tions very  similar  to  (1-3),  but  that  the  latter  cannot; 
the  same  is  true  about  these  climates  as  limiting  steady 
states  for  time -dependent  numerical  solutions  of  such 
equations . 

We  further  showed  how  changes  in  an  important  parameter, 
the  average  intensity  of  the  solar  radiation,  influence  the 
steady-state   solutions  of  the  model.   The  dependence  on  this 
parameter  of  all  steady  states  was  shown  to  be  gradual  and 
smooth  for  increases  of  up  to  5%  and  decreases  of  up  to 
about  2%.   However  for  a  critical  value  of  the  parameter, 
equal  to  98.15%  of  its  present  value,  the  glacial  and  inter- 
glacial climates  coallesced  and  they  disappeared  entirely  for 
smaller  values  of  the  parameter,  leaving  the  ice-covered  earth 
as  the  only  possible  stable,  steady  climate  of  the  model. 

This  result  is  important,  as  it  stresses  the  difference 
between  the  stability  of  a  steady  state  with  respect  to  the 
time  evolution  of  a  physical  system  governed  by  a  given, 
fixed  equation,  and  the  stability  of  a  steady  state  with 
respect  to  changes  in  a  parameter,  which  determines  the 
behavior  of  the  system.   For  definiteness ,  let  us  call  the 
former  internal  stability  and  the  latter  external  stability; 
we  have  shown  that  the  "deep  freeze"  is  stable  for  our  system 
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both  internally  and  externally,  that  the  glacial  is  unstable 
in  both  senses,  and  that  the  interglacial,  or  "present 
climate,"  is  internally  stable,  but  externally  unstable. 

The  limitations  of  equations  (1-3)  as  a  model  for  the 
description  of  the  long-term  behavior  of  the  atmosphere- 
ocean-cryosphere  system,  and  of  energy-balance  models 
in  general,  have  been  discussed  extensively.   Because  of 
these  limitations,  we  believe  that  the  results   above 
should  not  be  taken  at  face  value  as  statements  about  the 
climate  of  our  planet.   These  results,  however,  seem  to 
clarify  the  physical  content  and  mathematical  properties 
of  such  models.   Also,  the  methods  used  here  could  be 
helpful  in  investigating  other  models,  which  will  include 
more  elaborate  and  reliable  parameterizations  of  the 
physical  phenomena  governing  climate.   We  further  hope  that 
insight  gained  into  the  behavior  of  solutions  of  a  certain 
type  of  model  will  advance  the  formulation  of  other  models, 
and  that  these  will  come  closer  to  explaining  past  changes 
in  climate  and  predicting  future  changes. 
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Table  Captions 

Table  1.   Empirical  functions  appearing  in  equation  (4) . 
The  functions  Q,  b,  z   are  based  on  data  in  Tables  1 
and  2  of  Sellers  (1973).   The  functions   C,  k,  ,  k 
are  based  on  data  provided  by  Dr.  T.  Gal-Chen  (1974, 
personal  communication) ,  and  used  in  the  (SV)  model 
of  Schneider  and  Gal-Chen  (1973). 

Table  2.   Influence  of  different  modifications  in  the  model's 
equation  (4)   on  the  number  of  steady-state  solutions 
and  the  numerical  values  of  these  solutions.   The 
existing  solutions  are  identified  by  the  temperature 
at  the  pole,   u.  (0)  ,   j  =  1,2,3.   In  case  a  solution  is 
missing,  this  is  indicated  by  (-)  in  the  corresponding 
row-and-column  location.   S  stands  for  the  coefficients 
being  fitted  by  cubic  splines,  B  for  Bernstein  poly- 
nomials.  A  downward  arrow  (4-)  to  the  right  of  a  comment 
indicates  that  the  equation  used  in  the  numerical 
experiments  reported  in  all  subsequent  rows  had  the 
feature  pointed  out  in  that  comment.   Otherwise  comments 
refer  only  to  the  row  in  which  they  appear.   A  left 
arrow,   x  -*-  y ,  means  that  the  quantity  x  was  replaced 
in  the  equation  by  the  quantity  y.   The  entries  given 
with  less  than  five  decimal  digits  resulted  from  computa- 
tions with  lower  precision  than  indicated  in  the  text. 
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Table    3.      Dependence    of   the    steady   states    u^ (x; y) ,Up (x; p )    on    the 
normalized    average    intensity   of    the   solar   radiation,    y . 
The    columns    give      u(0),    the    temperature    at    the   pole, 
S   u,    the    average    temperature,    and      Au  =   u(l)    -    u(0) , 
the   pole-to-equator   temperature    difference    for   u,    and 
for   U2    respectively.      The    values    for      y      =    0.98151822 
correspond   to   the   bifurcating   steady   state      u   =   u    (x)  . 
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U2(0) 

Uj^(O) 

Comments 

168.94 

223.97 

247. 74 

S,  Sellers  a,  full  eq. ,  k  =  k(x,u) 

169.44 

227.76 

247.55 

S,   "           "       k  =  k(x,u) 

—             _ 

259.26 

S,   "        ,  Rq=  0.61au-^[u-4(u-u)  ]  , 

k  =  k (x,u) 

169.0 

222.6 

238.5 

S,   "                   k  =  k^  (x)  ^ 

170.0 

229.75 

238.04 

B,   " 

- 

229.2 

238.04 

B,   "                    a  /  0.85 

168.0 

— 

— 

S4-,  "         m  =  0  (no  greenhouse 
a<  0.85  i                                                effect) 

170.0 

222.0 

245.0 

Sellers  a,     c„=  0  i,  m  =  0.5  i 

170.0 

222.0 

- 

Faegre  a  4-,    a  /  0.2  5 

170.0 

222.0 

255.0 

a  ^  0.25  i 

169.0 

2  32.0 

265.0 

b(x)  ^  £^b(x)  =  2.85881  + 

- 

- 

- 

c(x,u)  -<-  0   (no  infrared  radiation) 

160.0 

- 

- 

V  cot  — ^  -^  0    (no  singularity) 
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250.0 

280.0 

k(x,u)  -r    £^k(x,u(x))  =  Kq  =  2.2xlo"^  i 

165.0 

- 

- 

V  cot  -^  ^  0 

- 

- 

- 

c  (x,u)  ■*-  0 

- 

247.36 

- 

Rq  =  a^c(x,u(x))ou'^(x)  =  5.63x10"^ 

193.27 

233.95 

277.79 

Rq  =  0.61  au^    (£^c(x,u(x))  =  0.61)  4 

190.0 

232.0 

280.0 

IT      .TTX     V 

V  2  cot   2  -  X 

197.0 

249.0 

295.0 

Q(x)  ^  £^Q(x)  =  8.333x10'^ 

190.0 

238.0 

275.0 

b(x)  ^  £^b(x)  +  Cj^C2£^z(x)  =  2.87334 

192.88 

232.08 

276.06 

k(x,u)  ^  K^  =  1.96x10"^                , 

Table    2 
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Figure  Captions 

Figure  1 

Comparison  of  curve  fitting  by  (i)  Bernstein  poly- 
nomial approximation,  indicated  by  a  dash-dot  line,  and  by 
(ii)  cubic  spline  interpolation,  indicated  by  a  solid  line. 
Bernstein  polynomials  are  not  interpolatory  and  they  are 
variation  diminishing,  i.e.,  they  have  the  property  of 
smoothing  out  the  data;  this  results  in  a  rather  poor 
approximation.   Cubic  splines  are  not  variation  diminishing 
and  they  are  very  good  approximants . 
Figure  la 

For  a  very  smooth  function,  like  u(x),  the  two  approxi- 
mation procedures  yield  curves  very  close  to  each  other. 
Figure  lb 

For  a  function  of  large  total  variation,  like  k(x,u(x)), 
the  two  procedures  yield  curves  which  can  differ  pointwise  by 
as  much  as  50  percent  of  the  average  value. 
Figure  2 

Numerically  obtained  values  of  u  (l;u„)   as  a  function 
of  Uq  =  u  (0)  . 
Figure  2a 

Comparison   of  the  results  for  equations  (7),  in 
which  k  =  k(x,u),  with  those  for  equations  (V),  in  which 
k  =  k(x,u(x));  the  results  for  (7)  are  indicated  by  a  solid 
line,  those  for  (7')  by  a  dash-dot  line. 
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Figure  2b 

Results  of  (V)  for  -1330K  '^  u{0)  <    300  K.   Notice  that 
as   u„  =  u(0)   tends  towards  the  ends  of  the  interval, 
u  (l;u^)  ^  +°°.   The  solution   u.(x),  corresponding  to  the 
negative  root  of  this  curve,   u^  =  -186  K,  does  not  have  a 
physical  significance. 
Figure  3 

Comparison  of  the  solutions  of  (4) ,  indicated  by  a  solid 
line,  with  those  of  (4'),  indicated  by  a  dash-dot  line.  The 
circles  indicate  mesh  data  for  u  =  u(x). 
Figure  3a 

Values  of  the  solutions   u . (x) ,  j  =  1,2,3,  for  (4)  and 

for  (4').   The  respective  values  for  j  =  1,3  are  practically 

indistinguishable,  whereas  for  j  =  2   a  slight  difference 

exists  between  the  solution   of  (4)  and  that  of  (4'). 

Figure  3b 

Values  of  the  derivatives  :^—   u.(x),  j  =  1,2,3,  for  (4) 

dx  J 

and  for  (4').   The  differences  are  larger  than  in  the  function 
values  themselves. 
Figure  4 

The  first  eigenf unctions ,   w.   (x)  ,  j  =  1, 2 , 3  ,  of 
L.  =  -  |-   N(u)  I 

1  dU  ' U=U . 

■'  D 

Figure  5 

The  function  u  (1;Uq),  obtained  by  integrating  (7') 
numerically  with  y  =  y   and  with  different  values  of  Uq , 


c 


u„  =  u(0)  >  120  K. 
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Figure  6 

Dependence  of  the  solutions  u.{x;y),  j  =  1,2,  on  the 
parameter  y.   The  two  plots  for   u.(0),   ^a^-;  ^^s  very  simi- 
lar; the  plot  for   Au  .  =  u.{l)-u.(0)   is  rather  different, 
although  it  exhibits  the  same  behavior  in  the  neighborhood 
of  the  critical  point  c.   The  circles  indicate  the  values 
actually  computed,  for  y  =  y   ,  0.982,  0.985,  0.99,  1.00, 
1.01,  1.02,  1.03,  1.04,  1.05.   The  letter  c  distinguishes 
the  values  of  the  plotted  quantities  u  .  (0)  ,  ^^u.  ,  Au. 
corresponding  to  the  bifurcating  solution  u  =  u  (x) . 

Figure  7 

The  bifurcating  solution  u  =  u  (x) .   Notice  that  the 
^  c 

ice  line,  which  corresponds  to  u  =  273  K,  is  at  about  45°  lat. 
for  this  solution,  i.e.,  an  ice  cover  extending  beyond  this 
latitude  would  eventually  cover  the  entire  earth. 
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